/*
 *  Klasa Model
 *
 *  Zawiera dane i metody dostepowe modelu.
 *
 *
 */
#include "model.h"

/*
 * Model::Model(int n) {
 * dlugoscLancucha = n;
 * W_[POPRZ].resize(n);
 * W_[NAST].resize(n);
 * P_.resize(n);
 * for(int i=0; i<n; ++i) {
 * W_[POPRZ][i] = W_[NAST][i] = 0.5 * multiplier;
 * P_[i] = 0.;
 * }
 * P_[0] = 1.;
 * }
 */

Model::Model() {}

Model::Model(string c) {
    lancuch=c;
    dlugoscLancucha=c.size();
}

int Model::getDlugoscLancucha() {
    return dlugoscLancucha;
}

string Model::getLancuch() {
    return lancuch;
}

double Model::getP(int i) {
    return P[i];
}

double Model::getW(int i, int j, double t, int zaburzenie) {
    //if(w[i][j]<1.e-6) cout<<"problem: "<<i<<" "<<j<<"\n";
    return 2*Pi*ws[i][j]*ws[i][j];
}

void Model::wczytajLancuch(istream iss) {
    iss>>lancuch;
    cout<<lancuch;
    dlugoscLancucha = lancuch.size();
}

void Model::przeliczSprzezenia(double t, int zaburzenie) {
    
    oscylacje = vector<vector<double> >(iloscBaz, vector<double>(iloscBaz, 1.));
    if(zaburzenie) {
        oscylacje[0][1] = calkaNakrywania(t, 0, 1, zaburzenie);
        for(int i=1; i<iloscBaz-1; ++i) {
            oscylacje[i][i-1] = calkaNakrywania(t, i, i-1, zaburzenie);
            oscylacje[i][i+1] = calkaNakrywania(t, i, i+1, zaburzenie);
        }
        oscylacje[iloscBaz-1][iloscBaz-2] = calkaNakrywania(t, iloscBaz-1, iloscBaz-2, zaburzenie);
    }
        
    ws = vector< vector<double> >(iloscBaz, vector<double>(iloscBaz));
    
    /* "zwykly" transport */
    ws[0][1] = w[0][1] * oscylacje[0][1]; 
    for(int i=1; i<iloscBaz-1; ++i) {
        ws[i][i-1] = w[i][i-1] * oscylacje[i][i-1];
        ws[i][i+1] = w[i][i+1] * oscylacje[i][i+1];
    }
    ws[iloscBaz-1][iloscBaz-2] = w[iloscBaz-1][iloscBaz-2] * oscylacje[iloscBaz-1][iloscBaz-2];
       
    
    /* tunelowanie */
    /*
    for(int i=0; i<iloscBaz-1; ++i)
        if(lancuch[i]=='G' && lancuch[i+1]!='G') {
        for(int j=i+2; j<iloscBaz; ++j) {
            if(lancuch[j]=='G') {
                
                /* w prawo 
                ws[i][j] = ws[i][i+1] * ws[j-1][j] / dE[i][i+1];
                for(int k=i+1; k<j-1; ++k)
                    ws[i][j] *= ws[k][k+1] / dE[i][k];
                
                /* w lewo 
                ws[j][i] = ws[j][j-1] * ws[i+1][i] / dE[j][i-1];
                for(int k=j-1; k>i+1; --k)
                    ws[j][i] *= ws[k][k-1] / dE[j][k];
                
                i = j-1;
                break;
            }
        }
        }
    */
}

double Model::calkaNakrywania(double t, int i, int j, int zaburzenie) {
    double c = 0.;
    if(zaburzenie == 2) { // WATER
        for(int k=0; k<3; ++k)
            c = c + ampW[k] * ( cos(freqW[k]*2.*Pi*t+fazaPocz[i]) - cos(freqW[k]*2.*Pi*t+fazaPocz[j]) );
    }
    else if(zaburzenie == 1) { // VACUUM
        for(int k=0; k<2; ++k)
            c = c + ampV[k] * ( cos(freqV[k]*2.*Pi*t+fazaPocz[i]) - cos(freqV[k]*2.*Pi*t+fazaPocz[j]) );
    }
    
    return exp(-1./(4.*sigmaSq)*c*(c+2.*DST));
}

void Model::inicjuj() {
    
    iloscBaz = dlugoscLancucha;
    
    w = vector< vector<double> >(iloscBaz, vector<double>(iloscBaz));
    dE = vector< vector<double> >(iloscBaz, vector<double>(iloscBaz));
    P = vector<double>(iloscBaz, 0.);
    fazaPocz = vector<double>(iloscBaz, 0.);
    
    P[0]=1.;
    
    vector<int> temp(iloscBaz);
    for(int i=0; i<iloscBaz; ++i) {
        switch(lancuch[i]) {
            case 'G': temp[i]=G; break;
            case 'C': temp[i]=C; break;
            case 'A': temp[i]=A; break;
            case 'T': temp[i]=T; break;
        }
    }
    
//    iTot=iN=0.;
    
    int s;
    for(int i=0; i<iloscBaz; ++i) {
        for(int j=0; j<iloscBaz; ++j) {
            w[i][j] = sprz[temp[i]][temp[j]] * multiplier;
            dE[i][j] = E[temp[i]] - E[temp[j]];
            dE[i][j] *= (dE[i][j]>0. ? 1. : -1.);
        }
    }
//    w[1][0]=0.;
//    w[iloscBaz-1][iloscBaz-2]=0.;
    
    for(int i=0; i<iloscBaz; ++i) {
        for(int j=0; j<iloscBaz; ++j)
            cout<<w[i][j]<<" ";
        cout<<endl;
    }
        
    srand( time(NULL) );
    for(int i=0; i<iloscBaz; ++i) {
        fazaPocz[i]=(double)rand()/RAND_MAX * 2.*Pi;
    }
}
